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7 Background of the Invention ^cS^o 

8 The present invention is directed to a method, program and 

9 apparatus for removing stray-flux effeGts in an image. * »-< 

10 All imaging systems are victimized by phenomena that misdirect a 

11 small portion of the entering flux to undesired locations in the image. The 

12 term "imaging systems" as used herein includes, for example, optical 

13 systems, x-ray systems, and computerized tomography systems. Other 

14 systems that do not employ the directional passage of flux, suchjas 

15 magnetic resonance imaging systems, also fall within the. term "imaging * 

16 systems"; Depending on the type of imaging system, the "image" may 

17 »- have a single dimension (a linear image), two dimensions (a. planar 

18 image), three dimensions (a volume image), or more dimensions (a hyper 

1 9 image): In the most general sense the misdirected flux may be termed 

20 "stray-flux", although in the context of different systems, it is known by 

21 I different terms. 

22 Optical imaging systems are victimized by phenomena that 

23 misdirect a small portion of the enteringlight flux to undesired: locations in 

24 the image plane. Among other possible causes, these phenomena 

25 include: 1) Fresnel reflection from optical-element surfaces, 2) diffraction 

26 at aperture edges, 3) scattering from air bubbles in transparent glass or 

27 plastic lens elements, 4) scattering from surface imperfections on lens 

28 elements, and 5) scattering from dust or other particles. The misdirected 

29 flux, which is called by the terms "stray light", "lens flare"* "veiling glare", 

30 and by other names, degrades both the contrast and the photometric 

3 1 accuracy of the image. For example, in photography, back-lighted scenes 

32 such as portraits that contain a darker foreground object, suffer from poor 

33 contrast and reduced detail of the foreground object. 
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1 Perhaps^ss known and appreciated is the eftt that "stray light" 

2 has on color accuracy of an image. P. A. Jansson and R. P. Breault, in 

3 "Correcting Color Measurement Error Caused by Stray Light in Image 

4 Scanners" published in Proceedings of The Sixth Color Imaging 

5 Conference: Color Science Systems and Applications, 1 998, Society of 

6 Image Science and Technology, referred to a traditional simple 

7 manipulation of offset and gain, or contrast, that can lessen subjective 

8 objections to stray-light contaminated images. These manipulations, 

9 however, do not correct the underlying flaw in the image and do, in fact, 

10 introduce additional error. 

11 U.S. Patent 5,153,926 (Jansson et al.), assigned to the assignee of 

12 the present invention, describes various-embodiments of a method to 

13 remove the stray-light effect from images. This method demands v 

14 significant amounts of computation which might inhibit application within 

15 an image acquisition apparatus, such as a digital camera. , l-^->: 

16 The method of the present invention, however, may be 

17 implemented in relatively simple apparatus, such as one or more digital- 
is logic integrated circuits within a digital camera, and can quickly and 

19 -inexpensively correct images that have been degraded by stray fight, 

20 Magnitude of the Stray Flux Effect The importance of the; stray-flux 

21 effect in general and the stray-light effect in particular may appreciated 

22 from the following discussion of the magnitude of the effect of stray light 

23 on an optical imaging system. In an image arrangement 10 as shown in 

24 Figure 1 , a transparent mask 2 having an opaque spot 3 may be. placed on 

25 a light-box 1 having a source of uniform diffuse illumination 4. A lens 5 is 

26 then .used to image this back-illuminated spot 3 and. its surround onto a 

27 target 6 comprised of a charge-couple-device (CCD) array of 

28 phbtodetector elements 6E, such as that now found in a digital camera. 

29 This arrangement creates an opportunity for light from the surrounding 

30 diffuse source 4 to stray Into the center of the spot 7 imaged on theTarget 

31 6. An opaque spot 3 having a lateral dimension of 1 2 percent of the 

32 surrounding field of view thus covers about 1 .44% of the image area. An 

33 example of an implementation of the image arrangement 1 0 has shown 

34 that the center of such an opaque spot 7 image on the target plane 6, that 
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1 should received light flux, actualiy receives light flWhat corresponds to 

2 about 2% of the adjacent surround. An idealized densitometer should 

3 indicate infinite (or indeterminate) optical density at the center of the spot 

4 in the absence of the stray flux. However, the 2% stray-contaminated data 

5 obtained from this imaging system results in an erroneous optical density 

6 value of D = log (1.0/ 0,02) « 1.7: 

7 Area-imaging systems collect light from all points in the image in 

8 parallel, in contrast to scanning imaging systems that collect light ,f rom 

9 each point in the image serially. From the above example, it may be 

10 appreciated that area-imaging optical systems have limited use when 

1 1 applied to the measu rement of optical densities. Best practice in optical 

12 density measurement calls for physically restricting the number of paths 

13 accessible to stray flux by using apertures, masks, baffles, special 

14 coating's and other techniques known to optical designers?,; ■<.. 

15 Scanning methods that image a scene in a point-byrpoint manner' 

16 can overcome the inaccuracy of single-exposure entire-image parallel 

17 detection. The serial acquisition of the image by mechanical scanning 

18 methods, however, is impractical for many applications, especially 

19 photography. The need to place an aperture or mask close to the object ' 

20 plane imposes another unacceptable physical limitation for photography. 

2 1 The greatly increased exposure time owing to loss of, parallel detection 

22 further imposes an unacceptable time limitation for photography. A 

23 photography camera requires a degree of separation between the object 

24 and the optics/detector system and requires parallel image detection to be 

25 practical.- -,• ••«>», >. • >....... ^ , • 

26 Color Accuracy in Photography and Imaging Colorimetrv P. A. 

27 Jansson' and R. P. Breault, in the earlier-cited publication consider the 

28 image arrangement of Figure 1, in a situation where a colored image spot 

29 is surrounded by a region having a different color, to determine how the 

30 color at the center of such a spot is affected by stray light from its 

31 surround. The perceptually uniform CIE color space denoted by the 

32 L*a*£>*coordinate system is employed." In this system a difference value 

33 of about 1 .0 is considered a just-noticeable-difference to a typical human 
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1 observer. Ablu^pot having 5% contamination of g^ren from a surround 

2 is shown to exhibit a color error of AEz.* a *b* = 8.94. Here the blue was 

3 specified to have the L*a*b*coordinates (10, 10, 30) and the green 

4 surround had L*a*b* coordinates (10, 30, 10). This shift of nearly nine 

5 just-noticeable-difference units, however, was small compared to the 

6 result computed when more saturated values of blue, having L*a*b* 

7 coordinates (6.5, 2.5, 40), and green, having L*a*b* coordinates (6.5, 40, 

8 2,5), were used, even with a much lower contamination level of 2%. In 

9 this case, the color error was gross, giving A£/.* a */>* = 1 5.07. Similarly, in a 

10 realistic black and white image, they determined that a 4% contamination 

11 of a black spot having 5% reflectance resulted in an error of AElw ■= 

12 9.09. ■ - .v • .. ■;, - - 

13 : ; Earlier-cited U.S. Patent 5,1 53,926 may be applied to the removal 

14 of stray-light effects from scanned images acquired by densitometers, 
15r_; cameras, microscopes, telescopes and other optical imaging systems. 

16 ! The method of this: patent contemplates computations in .both the Fourier 

17 r? and spatial domains. They require the concept of a full point^spread 

18 function (PSF) that characterizes the redistribution of light flux in the image 

19 . plane occasioned; by the presence, of scattering phenomena and'other ; 
20. stray-light-pfoducing effects in the optical system. Gjiaraeterizatpirpf i& 

21 stray-light by a point-spread function is not typically done because it 

22 describes the spreading of flux over large distances of image: coordinate, : 

23 This contrasts with common usage in which the point-spread function 

24 incorporates primarily diffraction and geometrical-optics aberrations. Only 

25 the central core of a camera's full point-spread function would correspond 

26 to the more traditional-point-spread function. 

27 . However, the following description considers'the point-spread 

28 function to provide for replication of an object with added stray-light effects 

29 only. Moreover, the stray light distribution is typically shift-variant. A 

30 single convolution is not usually adequate to describe stray-Hight 

31 contamination over a complete image. Accordingly, imaging in a digital 

32 . camera is described by specifying all of the image flux 7 (x, y) received by- 

33 a rectangular planar image detector having surface spatial coordinates x 

34 and y: 
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1 i(x,tf = tis(m;x,y')o(x,y)dxdy. ♦ (1) 

2 In this equation, x and y are the two-dimensional spatial coordinates of a 

3 planar object o {x, y) that is imaged onto the detector, and s (x, y; x, y) 

4 is the point-spread function that describes replication of the object with 

5 only superimposed stray-light effects, but no additional corruption. The 

6 object distribution o(x', y) according to the present definition therefore 

7 may exhibit the other effects typical of an imperfect imaging system, such 

8 as degraded resolution owing to diffraction and geometrical aberrations. 

9 The present method also assumes that all of the flux arriving at 

10 image /(x, y) originates within, the bounds of the object that map directly to 

1 1 the rectangular bounds of the detector, so the limits of the integral are 

12 specified by these object bounds. 

13 Jansson and Breault, in the previously cited publication,, had shown 

14 that,; because the character of point spread functions, such as ■;, ■*■ 

15 . :,;$;<*, y, . xr, y% necessarily vary with the wavelength of light, both the 

1 6 imaging-system description and stray-light correction methods need to 

17 consider this variation in order to rigorously correct color images. The 

1 8 present method, in order to simplify both camera description and image 

19 correction, considers the CIE standard-observer j> and ~z functions as 

20 each describing separate monochromatic responses with wavelengths 

21 corresponding to the peaks of these functions, respectively. Following the 

22 same reasoning strictly, the x function would require summing responses 

23 at the two wavelengths that correspond to its two peaks. However, color 

24 cameras typically employ detectors having a single-peak wavelength 

25 response. For clarity, notation of wavelength dependencies may be 

26 suppressed, while noting that required computations can be replicated at 

27 numerous multiple wavelengths for completeness, or at detector/receptor 

28 response peaks as a practical approximation. 

29 Rigorous correction of image / (x, y) for stray-light effects requires 

30 solving Equation (1 ) above for object spatial distribution .o (x, y). 

31 Problems similar to this.one, but in which s (x, y, x',y) is shift invariant, so 

32 that s (x, y, x, y) can be simplified as s (x - x, y- y), call, for 



5 



1 deconvolution Clause Equation (t) would then con^to a two-dimensional 

2 convolution, i.e., 

3 i (x,y) = JJ s (x - x\ y- y) o (*', y) dx'dy. (2) 

4 It is typically difficult to solve such problems involving real-World data 

5 inputs because they are usually ill-posed. By this it is meant that they 

6 either have no solution, have numerous solutions but offer no practical 

7 means to distinguish amongst them, or- have a solution that varies radically 

8 or discontinuous^ with small changes in data / (x, y). It might be expected 

9 that the present shift-variant problem would be similar. However, this 

10 turns oiit not to be the case, for reasons'that become apparent when the 

11 discrete version of Equation (1) is considered. ■ ' 

12 Consider then the matrix equation 

13 i = so (3) 

14 that captures the sampled version of the problem, the column matrices i ! 

15 and 6 corresponding to / (x, y) and o (x', y), respectively, and the square 

16 matrix s to s (x, y, x,y). The double-variable X and y spatial 

17 dependencies of / (x, y) and o (x', /)■ collapse neatly to the single variable 

18 arguments required by the matrix description via lexicographic ordering of 

19 object and image planes within their bounds. See E. J. Borowski and J. 

20 M. Borwein, The Harper Collins Dictionary of Mathematics, Harper . 

21 Perennial, New York, 1991 , p. 339, Accordingly, the dimensions of 

22 matrices land o are given by m x n, ? where m is the number of samples 

23 across the width of object o (x, y) and image /'(x, y), and n is the number 

24 of samples along their respective heights. Matrix s then has dimensions N 

25 x N, where N = m x n. 

26 One method of solving for vector o calls for iterations of the equation 

27 6 ( * fl) =6 ( * ) + (i-s6 ( * ) ), (4) 

28 where k is the sequence number of the approximate solution, © w is the lP 

29 estimate of the uncohtaminated object and i is typically taken as 6 (0) . For 

30 the shift-invariant case involving discrete convolution this' is called van ' 

31 Cittert's method. It has been identified as a special case of the Jacobi 

32 method of solving a set of linear equations, which in turn is applicable to 

33 the more general shift-invariant case applicable here. See P. A. Jansson, 
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1 "Traditional LJPar Deconvolution Methods," In DecWvolution.of Images 

2 and Spectra, P. A. Jansson, ed., Academic Press, New York, 1 996, p. 85- 

3 89. 

4 In order to show that the solution to Equation (3) is not ill-posed, 

5 consider some specific features of s. In particular, let the PSF of a perfect 

6 camera, devoid of stray-light contamination, be denoted by s P = I, where I 

7 is the identity matrix. The stray-light component may be characterized by a 

8 matrix So, the sampled version of. an equivalent continuous function 

9 s D (x, y, x', y) so that ... 

10 s = Sp+s D . (5) 

11 Because the present concern is primarily with stray light, either of 

12 three assumptions can be made: (1 ) that imaging-aperture diffraction and 

13 geometrical aberrations are negligibly small; (2) that aberrations are 

14 spatially localized well enough to be contained within the diagonal 

15 elements of the matrix sp, or (3) that such effects are already incorporated 

16 in object.,distribution described by vector o. Any of the three assumptions 

17 leads to the same mathematical solution. However, they differ in their 

18 influence on how the solution must be interpreted and used. In the case of 

19 the thjrd assumption, even a perfect solution will yield an object estimate 

20 6 bearing, the effects of an otherwise imperfect imaging system. /. 

21 Accordingly, the present method corrects for neither airy-disk-like 

22 diffraction owjng to a. finite lens aperture nor for lens aberrations. 

23 Note that in the foregoing, it has also been assumed that the scale 

24 of the object and image is such that a value of unity results along the 

25 diagonal. This assumption is necessary for the results that follow. 

26 Furthermore,.the individual sums along all rows and columns necessarily 

27 . exceed unity owing to the presence of stray light: 

28 I ,[s],> > 1 for all /, and I y [s] (> > 1 , for i. (5 A) 

29 Consider the first iteration solution, obtained by substituting the starting ,. 

30 estimate 6 (0) = i into Equation (4); 

31 6 (1) = i+.(i-si) (6) 

32 = (spp + s D o) + [ (S/.0 + s 0 o) - s P s P o - 2 SpS 0 p - s D s 0 o] ... , 

33 Noting that s p = 1- terms may be collected, thereby yielding 
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1 6 (l) = o-s D s D l^ ^ (7) 

2 This is an especially interesting result, a serendipitous consequence of the 

3 manner in which s P is normalized. The first iteration object estimate is the 

4 object itself, in error by merely the quantity - Sds d o . This quantity is 

5 small because the stray-light portion of a digital-camera image is small 

6 compared to the properly imaged portion so that 

7 lyfeiy « Z/[S P ],y. (8) 

8 This result obtains unless, of course, the camera's lens is grossly 

9 contaminated. 

10 As seen above, a useful approximation to the uncontaminated 

1 1 object haS been obtained by only one cycle of an otherwise iterative 

12 computation, this cycle involving a single matrix product: 

13 6 ;(1 > = 2i-si 

14 ^ i-sbi. ^ • ' : " (9) 

15 In-order to correct the image for stray-light contamination to good 

16 approximation, subtraction from the image of a term i comprising the 

17 convolution of the image with the point-spread function stray-light 

18 component may be performed. Jansson and Breault, cited above, used a 

19 similar method to correct one-dimensional images for simulated stray light, 

20 yielding excellent results in the case Where the PSF was characterized by 

21 a hyperbolically decreasing shift-invariant stray-light contribution 

22 superimposed on a "perfect-imagihg"naifow central Core. 

23 Even though the aboVe simplification has greatly reduced the 

24 potential computational load, it's direct extension to two-dimensional 

25 images nevertheless requires a large amount of computation if the sums 

26 are computed directly in the spatial domain.' Specifically/an 'N by N image 

27 would require A/ 4 multiply-and-add pairs of operations, or 2/V 4 operations 

28 altogether. 

29 If the point-spread function were shift invariant, the computation 

30 would involve a convolution that could be accomplished in the Fourier 

31 domain. If it is assumed that N is power of 2, a one-dimensional FFT of a 

32 single image line would require 5 N 2 log2 N operations. See P. A. Jansson, 

33 "Convolution and Related Concepts," in Deconvolution of Images and 
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1 Spectra, P. A.Wnsson, e<±, Academic Press, 1996^25. The 2-D 

2 transform then would require 2N times that number, or 1 0 /V s log2 N 

3 operations. 

4 One transform would be required to convert i to the transform 

5 domain, where it would be multiplied by-the stray-light transf etf unction 

6 (21 - s). This transfer function would already be in place, having been 

7 previously transformed from spatial-image space and available for multiple 

8 uses. A second transform would take the product back to the spatial 

9 domain. Thus the.convOlution involved in correcting a single image would 

10 cost two transform operations plus 6 /V s multiplications in the transform 

11 domain, or (6 + 20 log 2 N) operations altogether. 

12 However, it is already known that s D (x, y, x', y) is not shift-invariant. The 

13 correction can still be carried out via convolution by breaking the image ; 

14 into numerous small patches Over which s(x, y x, y') does not van/. 

15 significantly, so-called isoplanatic patches. Note that, although the*image 

16 data require transformation only once, convolution-theorem multiplication 

17 and inverse transformation must be carried out for all the patches. If it is 

18 • assumed the patches are M x M pixel squares, the total number of 

19 operations needed becomes > . k 

20 N 2 {I0i 1 * *(N/M) 2 ) log2 ^+ 6 (N/M) 2 y~ {N*/M 2 )i 101og 2 W^ 6 ) . (9A) 

21 This too imposes a large computational burden and requires computer 

22 code having increased complexity. Before introducing a preferred method 

23 of reducing this burden and complexity, it is appropriate to introduce a 

24 method to apply traditional multi-resolution approaches to efficiently 

25 compute the stray-light correction. u 

26 Image Pyramids It is well known that much computation can be 

27 saved in processing and analyzing images by operating On reduOed- 

28 sample-density images in which each pixel is the composite, for example, 

29 the average, of a group of pixels. See K. R. Castleman, Digital Image 

30 Processing, Prentice-Hall, 1996, p. 320. Once the rough features of the 

31 analysis have been obtained on the reduced resolution image, it Can be 

32 refined by reference to the fully sampled original. To this ehcl, it is even 

33 possible to compute a hierarchical sequence of images to accomplish a 
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1 given computaf^nal goal in stages! For example, colder a so-called 

2 pyramid of images in which each successively smaller (i.e., higher on the 

3 pyramid) image layer contains pixels corresponding to two-pixel by two- 

4 pixel square regions of the preceding image. Operating on just the first 

5 reduced-sample version of the original would result in a four-fold 

6 computational savings. Storing.a complete pyramid of such images 

7 requires' ortly /v* + (A// 2) 2 + (N/ A) 2 + . . . + 1 image pixel locations. 

8 The effect that such a sampling method has on the stray-light- 

9 correction computation may be considered. First a pyramid may be 

10 computed from the original image data /' (x, y). In correcting pixels where 

1 1 the So (x, y, x, y) varies rapidly, near the central core, data from the f ully- 

12 sampled image would be used. Farther from the core, so (x, y, x\ y") 

13 varies more slowly with distance. Pixels from a more coarsely sampled 

14 image in the pyramid may be employed to form the products needed, 

15 using the weights from this portion of s D (x, y, x\ /). Ideally an > 

16 appropriate level from the pyramid could be chosen such that the error 

17 dOntributions from all the samples-would be about comparable. 

18 Depending oh the specifics of the s D (x, y, x, y) and of the computational 

19 accuracy required, much computation can be saved via this method. The 

20 computational bookkeeping required to implement the method, however, 

21 could be substantial, 

22 The idea of varying the sample density leads to yet another 

23 approach. In the example above, the sample density varies by factors of 

24 two between layers of the pyramid. As the values of s D (x, y, x, y) are 

25 traversed, the pyramid's structure forces adjustment of sample density by 
26. a minimum factor of two. In consideration of the error contribution, the 

27 least advantageous samples can have an area-based density as much as 

28 two times greater or less than that which is ideal. This observation 

29 encourages examination of smaller sample-ratio increments. A pyramid of 

30 images having very small increments would contain a great many layers 

31 and occupy much storage memory, so this option might be dismissed as 
32. impractical, especially in the light of the more attractive alternative 

33 described in the following. 
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1 If a samPIng pattern is computed in which deHity varies locally in a 

2 smooth manner, such that it adapts to a uniform error-contribution 

3 requirement, the disadvantages of the pyramid method can be overcome. 

4 More significantly, in such a sampling pattern the sample density itself can 

5 be used to encode the weighting required by 

6 s D (x, y; x y '). Thus, explicit multiplication may be completely eliminated 

7 in computing the integral or matrix product. 

8 Selected Ordinate Concept Consider the. accuracy and precision 

9 requirements of computing an arbitrary one-dimensional weighted integral 

10 g = jw(x)f(x)dx, (10) 

1 1 where f (x) is a function to be integrated, w (x) supplies the required " 

12 weighting, and x is the independent variable of integration. The value of a 

13 two-dimensional integral product such as the one in Equation (2) may be 

14 obtained for any corresponding value of its argument by computing just 

15 such an integral, and the discrete product required for sfray-light'correction 

16 by Equation (9) offers a similar case. Equation (10) may be approximated 

17 by the discrete sum so that 

18 g-il/Nfc^fM, (11) 

19 where Axy= Ax 2 = . . . =Axat is the uniform spacing between N samples. 

20 In order to minimize the labor of computing this approximation, it is 

21 desirable to distribute the burden of precision equally amongst the terms 

22 summed, and equally to each multiply-and-add pair of operations; i.e., 

23 each term should contribute equally to the resulting error. If the variances 

24 of terms to-, f, are equal, a simple sUm will suffice to this end. However, 

25 this would be an unusual case, to say the very least. 

26 In a practical case, there might be no prior knowledge of f(x): 

27 Considering the precision of g in complete absence of such knowledge, a 

28 reasonable choice is that the // are equal. Under these circumstances 

29 terms wj f, do not contribute equally to the sum, so each multiply-and-add 

30 pair of operations contributes unequally to the precision of the result. 

31 Being confined to a set of evenly spaced samples f h there is no means to 

32 achieve the goal of error contributions that are uniform over a desired 

11 



1 domain. An alt^ate weighting of the terms would police a result that is 

2 not consistent with the definition of g. The clue to improving precision and 

3 computational efficiency is found by examining the assumption of equal 

4 sample spacing. 

5 In accordance with the foregoing, if the sample spacing Ax/ is 

6 adjusted so that each term contributes equally to the error, optimal use is 

7 made of the labor expended not only in sampling, but in computation as 

8 well: Specifically, let <the local sample spacings now vary inversely with 

9 weights w, such that spacing Ax, = 1 / w<, incorporates the effect of the 

10 weight. Here 1 / w, may be taken as the spacing between the (/- 1)* 

1 1 and the F samples, the P and (/ + 1 samples, or the spacing may be 

12 divided in some other manner in proximity of the sample. For example, 

13 the spacing between the P and (/+ 1)* samples may be taken as 

14 0.5 (( 1 / Wj) + ( 1 / Wj+i)). By any of these methods the weights no longer 

15 appear explicitly in the sum. Because computation of the integral no 

16 Jonger involves even a single multiplication, the number of arithmetic 

17 computations is immediately halved. Moreover, because the samples are 

18 placed only where they are the most needed, i.e., where the Wj has its 

19 largest values, fewer samples are required for a result having a given 

20 .j precision. In particular, Jn the ca$e of photon-counts or Poisson^detection 

21 noise, assumed uniform f (x), and a collection of all counts over the 

22 specified interval, each term on average contains the same number of ; 

23 counts, has the same variance, and thus contributes equally to the 

24 precision of g, which is the goal. 

25 Consider also the case in which f (x) is affected by additive noise. 

26 Here f (x) may be finely sampled with even sample spacing. The uneven 

27 spacings Ax, described in the foregoing may then be used to establish 

28 which of the fine samples are to be included in the sum, the other samples 

29 being ignored. In this case, as before, all the samples included in the sum 

30 have equal variance, though the information from the fine samples 

31 between those selected does not contribute to improving the precision of 

32 the sum. 

33 Before the advent of digital computers, spectroscopic colorimetry 

34 employed a technique of adjusting sample spacing inversely with the 
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1 desired weighwg to achieve efficient computation affhe weighted 

2 integrals. No multiplications are necessary with this selected-ordinate 

3 method. See D. B. Judd and G. Wyszecki, Color in Business, Science, 

4 and Industry, 3rd ed. J. Wiley and Sons, New York, 1975, p. 149. An 

5 operator records a spectral-reflectance curve on chart paper pre-printed 

6 with vertical lines, or ordinates, selected according to the weighting 

7 required by the standard-observer functions. The operator then merely 

8 reads the ordinate values at the selected locations and adds them to 

9 compute the required weighted integral. Though a primary motivation is to 

10 eliminate multiplication, far fewer samples are found Xo be necessary for 

11 accurate color determination. 
12 

13 Summary of the Present Invention 

14 The present invention is directed to a method, program and 

15 apparatus of correcting an image of an object, wherein the image 

16 comprises one or more dimensions and has a plurality of image elements 



17 and the image, or the data used to create the image, is corrupted by 

18 smearing or misdirected flux. The present invention removes the effect of 

19 the corruption on the image, the method comprising the steps of: 

20 a) sampling portions of the image of the object at unequal, intervals, 

21 the intervals being selected to effect an implicit weighting of the samples 

22 corresponding to a predetermined weighting function; 

23 b) summing the implicitly weighted samples to create an 

24 approximation to a weighted integral; 

25 c) repeating steps a) - b) to compute weighted integrals associated 

26 with each of the.plurality of image elements thereby to obtain a correction 

27 to the image. 

28 The present method further comprises the step of determining a 

29 weighting function that characterizes the corruption before performing step 

30 a). After obtaining the correction to the image the correction is applied to 

31 the image to remove the corruption and to create a corrected image. 

32 The method may be applied in an iterative manner and may further 

33 comprise the steps of: applying an additional correction term containing a 

34 sum of implicitly weighted samples to each element of the corrected image 
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1 to update the ejected image and then repeating th^tep until a 

2 predetermined stopping criteria has been satisfied. 

3 The method of the present invention may also be applied to 

4 situations where the number of dimensions of the image differs from the 

5 number of dimensions of the object. 

6 In a preferred embodiment, the method of the present invention 

7 uses sampling patterns in two dimensions analogous to the one- 

8 dimensional, hand performed, spectroscopic sampling to compute the 

9 weighted integrals necessary to apply methods such as those expressed 

10 by Equation (4) and Equation (9). 

11 ' ■ - 

12 Brief Description of the Figures 

13 The invention will be more fully understood from the following 

14 detailed description, taken in connection with the accompanying drawings, 

15 in which; < 

16 1 Figure 1 shows an image arrangement which illustrates the effect of 

17 stray light on an image; 

18 Figure 2 presents pseudo code showing the sample-coordinate 

19 selecting process; •' * 

20- Figure 3 presents pseudo. code showing the process of cqmputing a 

21 weighted integral using the selected ordinate selected coordinates; : 

22 Figure 4 shows an apparatus which may be used to implement the 

23 present method; 

24 Figure 5 depicts a highly simplified case of selected ordinate 

25 sampling of the image within an image; 

26 Figure 6 shows the image samples selected to correct an adjacent 

27 pixel in the same isoplanatic patch of the image of Figure 5. 

28 

29 Detailed Description of the invention 

30 Throughout the following description, similar reference characters 

3 1 refer to similar elements in all Figures of the drawings. 

32 Generating the Predetermined Weighting Function Before the 

33 selected ordinate image sampling of the present invention is performed, a 

34 step of determining a weighting function that characterizes the corruption 
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1 to the image ipjerformed. For applications to objeW and images having 

2 two-dimensions, the weighting function may be generated by any of the 

3 methods outlined in U.S. Patent 5,153,926 for producing the response 

4 coefficients required therein, or by any other method suited to this task. In 

5 U.S. Patent 5,1 53,926, the imaging apparatus is calibrated using a test 

6 object having known properties. This test object may be the preferred 

7 object comprising a plurality of flux sources, singly excited in some 

8 predetermined order, or another known test object, such as a high- 

9 contrast edge. In a publication authored by Berge Tatian "Method for 

10 Obtaining the Transfer Function from the Edge Response Function," ; 

1 1 Journal of the Optical Society of America, Volume 35; Number .8;. August 

12 1 965, the use of a high-contrast edge test object; is described.: It should be 

13 appreciated that the transfer function may be Identified with the Fourier 

14 transform of the desired weighting function. Also, it should be appreciated 

15 that other forms of test object may be found useful, or the response 

16 coefficients may be calculated based upon the design of the imaging 

17 apparatus, or the Weighting function may be selected arbitrarily. 

18 By analogy to the two-dimensional imaging case, weighting 

19 functions required by imaging systems having other dimensionality, or by 

20 imaging systems in which object and image have differing dimensionality s 

21 may be determined by use. of a suitable known test object. 

22 Generating the Selected-Ordinate Pattern In accordance with the 

23 present invention, 46 sample an image according to a Specified weighting 

24 function using selected ordinates, first sample locations are generated 

25 such that the local sample density varies in accordance with the weighting 

26 function. The resulting , sample pattern depends on how the "local sample 

27 density" is determined. In a first "local sample density" method, the 

28 average distance from a given weight location to the P closest neighboring 

29 samples is computed. The square of this distance then becomes a 

30 measure of the area covered by the samples. It is recognized that 

3 1 geometric distribution of the P samples influences the precise area 

32 covered, but for practical purposes this effect may be ignored. A first 

33 alternative "local sample density" method may explicitly build the precise 

34 area into the density estimate. For example, the number of samples within 
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1 a specific geomPh'cally defined neighborhood may bWounted. In a 

2 second alternative method the sample contributions to the density 

3 estimate could be weighted according their relative position within the 

4 neighborhood, thereby giving the neighborhood "fuzzy" boundaries. 

5 To apply the method of sample-pattern generation, the local sample 

6 density, according to whatever definition is used, at each pixel is 

7 compared to the weight desired for that pixel. If the density is smaller than 

8 that specified by the weight, or a fixed predetermined fraction of that 

9 weight, a sample is assigned to that pixel location. In practice, it has been 

10 foundaiseful to assign the first sample at the location of the weighting 

11 function peak, and to iterate (cycle) through all of the image coordinates 

12 several times. During the first cycle, an attempt is made to satisfy only a 

13 small fraction of the targeted weighting. Subsequent iterations (cycles) 

14 satisfy larger and larger fractions until the last pass, at which time the full 

15 weight sample densities are obtained. The result is aJull complement of 

16 the sample locations needed to achieve the desired weighting. The 

17 sample-coordinate selecting process is illustrated by the pseudo-code 

18 listing Of a program for implementing the method of the present invention 

19 shown in Figure 2. , 

20 In order to balance the requirements of computational load against 

21 sampling accuracy? the total number of samples may be constrained. 

22 Also, the result of implementing a computation, such as that of Equation 

23 (.1 1), might yield a value that is scaled smaller or larger than the desired 

24 value, yet otherwise exhibits the features that are desired with respect to 

25 weighting. In this circumstance, the results of all the pixel computations 

26 for a single isoplanatic patch may each be multiplied by a single constant 

27 scale factor, determined for that patch, thereby accommodating scaling 

28 imposed by a number-of-samples constraint, sample sparsity, and other 

29 phenomena that affect scaling. > ? 

30 In some circumstances it may be desirable to use two or more sets 

31 of samples, adding their contributions with different scale factor 

32 weightings. 
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1 The apflkach of adaptively varying the inter^between image 

2 sample coordinates in the present context bears a certain appearance of 

3 similarity to methods used to place half-tone dots for printing images. In 

4 the present case, it is desired to minimize the error in the image by optimal 

5 choice of sample coordinates. In the half-tone dot case it is desired to 

6 maximize the acceptability of appearance of a displayed image by optimal 

7 choice of the half-tone dot coordinates; In the present case of selecting 

8 image sample coordinates, the ideal function is related to the. point spread 

9 function s. In the half-tone dot case the ideal function is the image to be 
10 1 displayed and the dot coordinates are not used to sample another 

1 1 function, as in the present case. For a.discussion of half-toning; see D. L. 

12 Lau and G; R. Arce, Modem Digital Halftoning, Marcel Dekker, New York, 

13 .2001. ..... .... 

14 •.. .. . . 

15 : ;j- Isoplanatic Patches If theiorm of the point-spread function is such 

16 , that the convolution of Equation (2) is an adequate approximation to the 

17 more general form of integral product of Equation (1) over a limited region 

18 of an image, such a region is called an "isoplanatic patch". In this case a 

19 single set of weights may be employed over the patch by adding x and y 

20 offsets to the selected ordinate sample locations as they are used to 

21 sample the image. 

22 The general point-spread-function kernel may be re-formulated as a 

23 new function containing a set of x and / difference coordinates by defining 

24 s g (x - x\ x'; y- y', y) = s (x, y, x, y). , (12) 

25 An isoplanatic patch may then be defined as one over which 

26 s s (x->x'., x'; y - y\ y) varies slowly with x'and y 'even though it may 

27 vary rapidly in the difference coordinates x- x'and y- y'. Within that 

28 patch we then can have 

29 s b (x-.x\yry) = Sg{x-x,Xo,y-y',y 0 ), (13) 

30 where x 0 and yb may, for example, be patch center coordinates. Each 

31 patch has, of course, a different s b . 

32 , Accuracy of the Selected. Ordinate Estimate In order to roughly 

33 assess the number of samples required to correct an image for stray-light 

34 contamination, consider an image in which each pixel is contaminated with 
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1 "salt?and-pepp^bimodally distributed noise rangin^wer 5% of full scale, 

2 i.e., from -/ 0 / 40 to +/ 0 / 40. Specifically, the image pixels are all affected 

3 by error contributions having values of either -/ 0 / 40 or +/ 0 / 40. These 

4 contributions siffect each and every image pixel, and are present in equal 

5 number for each type of contribution, the type of contribution affecting a 

6 given pixel being randomly chosen.: The variance of this noise is given by 

7 Iq 1 1600. If a normal distribution of the noise error is assumed, the 

8 variance of the of /.-sample mean is given by 7 0 2 / 1 600L . Employing 64 

9 samples to estimate the stray-light contribution at each pixel implies that 

10 the estimate will have standard deviation (/ 0 2 7 1 02400 ) 1/2 « 0.0031 /o- If it 

1 1 is anticipated that stray-light contaminates each pixel at a level of at least 

12 several percent, the image would seem to be adequately sampled via the 

13 selected ordinate method to provide an effective estimate of stray4ight 

14 content and therefore should afford substantial improvement. 

15 * ; It should be noted that no computational method of correcting for * 

16 stray-light contamination based oh the image data alone can accurately 

17 correct for contamination by sources of flux falling outside the bounds of 

18 the image. This is caused by a lack of data required to determine the 

19 degree of that "out of bounds" contamination. Such "out of bounds" 

20 contamination may be estimated from data within the image. Thjs 

21 estimate may be used to approximately correct for such "out of bounds" 

22 contamination. Lenses may be designed with features that minimize, but 

23 do not totally eliminate, image contamination by stray light from such 

24 sources. 

25 Computational Advantage of Selected Ordinate Method Assume 

26 that it has been found sufficient to sample the image at L selected 

27 ordinate^ in order to achieve a selected-ordinate-based estimate having 

28 the desired statistical precision at each image pixel. Because only 

29 addition operations (no multiplications) are needed, approximately Llf 

30 operations will suffice to process the entire image regardless of the 

31 number of isoplanatic patches. Combining this result with the expression 

32 on Line (9A) above, it may be determined that the selected-ordinate 

33 method runs faster than the Fourier method by a factor of 
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1 (10log 2 N + ^fc 2 /M 2 L. • W (14) 

2 To see whether improvement is significant in a practical case, 

3 consider typical values for the variables in this expression. Consider a 

4 four megapixel image, having dimensions 2048 x 2048, broken into 64 

5 isoplanatic patches, each patch having dimensions 256 x 256. Previously 

6 it was found that 64 selected-ordinate samples were sufficient for stray- 

7 light correction to within 0:31% in a-representative case. Substituting the 

8 values N = 2048, M =. 256, and L = 64 into the expression in Equation (14) 

9 gives a 11 6-fold improvement over an FFT-based method — over two 

10 orders of magnitude reduction in computation time. 

11 Pyramidal Selected-ordinate Sampling Selected-ordinate 

12 sampling might be expected to work especially well for images containing . 

13 distributions of flux that vary slowly with coordinate and distributions over a 

14 - small dynamic range. However; images containing point sources, for 

15 example/ would be problematic because a given sampled pixel might 

16 pooifyrrepresent its neighbors. For these cases, the pyramid concept 

17 comes to the rescue. Where the density is low, the sample can be chosen 

18 from a high level in the pyramid, where each value already represents an 

19 average of neighboring values. Conversely, where the sample density is 

20 high, a lower level would be chosen. ■ 

21 Efficient Computation and Hardware Once a sampling pattern has 

22 been established, the locations s>f the samples for a given isoplanatic 

23 patch can easily be stored in a lookup table. Stepping through the lookup 

24 table generates the coordinate locations in sequence. Pixel values stored. 

25 at those locations are merely summed to produce a weighted sample of 

26 the image that requires at most only a single multiplication for scaling 

27 purposes. When applied to a convolution product over a particular 

28 isoplanatic patch, the coordinate pair of values from the lookup tables is a 

29 pair of offsets to be summed with the coordinates that locate that particular 

30 patch. The result of this pair of sums gives the sample location relative to 

31 the full image. A psuedo-code listing that applies the selected ordinate 

32 samples to calculate a convolution product over a particular isoplanatic 

33 patch in accordance with the present invention is shown in Figure 3. 
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1 ImplemeWng a convolution over a finite spatiarclomain involves 

2 consideration of the special treatment kernel samples must be given as 

3 they shift beyond the bounds of that domain. The effect is particularly 

4 pronounced when the kernel is broad. In the case of strayrlight, this 

5 kernel, So, extends over the entire domain. Because the image is to be 

6 divided into patches, however, special processing-is required only in the 

7 so-called "border region" of the image. The breadth of the border region 

8 is nominally half the patch breadth, assuming the patch's (xo.yo) pixel is 

9 centered in the patch. Minimizing this border's breadth is another good 

10 motivation to divide the image/into such patches. A preferred special 

1 1 treatment is the substitution of the closest in-bounds pixels for samples 

12 specified beyond the bounds, Of course another suitable approximation 
13. could be used. ■ ■ ■■>. ■ 

14 :• Note that throughout the present method, no Fourier transforms are 

15 required, arid no multiplications are needed in' the inner computational 

16 loops. Simple integer arithmetic sums are sufficient. The computation is 

17 sufficiently simple that design of special-purpose processor hardware 

18 specialized to this task is straightforward. The principal features of such a 

19 processor are outlined by the block diagram Figure d. • 

20 Figure 4 shows an apparatus generally indicated by Jhe reference < . 

21 character 1 00 for correcting art image of an object in accordance with the 

22 method s the present invention. The apparatus 100 comprises an image 

23 memory 1 10 having one or more partitions^ 1 0A, 1 1 0B for storing an input 

24 image and a corrected image; an arithmetic logic unit 120, comprising sub- 

25 units 120A, 120B, 120C, 120D, for performing arithmetic computations; a 

26 plurality of look up tables 1 30A, 1 30B, 1 30C, . . > 1 30N for storing selected- 

27 ordinate image coordinate values; a scale factor look up table 1 40, a 

28 controller 1 50; comprised of sub-units 1 50A, 1 50B, 1 50C, 1 50D, and an 

29 image address generator 1 60. The functional elements 1 1 0, 1 20; 1 30, 

30 140, 150, 160 may operate under a program in accordance with the 

3 1 present invention for addressing the image memory 110 and selecting 

32 lookup tables 130. The program may be stored on any suitable computer 

33 readable medium. Alternately a full hardware implementation of the 

34 functional elements may be employed. 
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1 In opeflfcn, first an image of the object is sflkd in the image 

2 memory 1 1 0A. The controller 1 50C, under the control of the program, 

3 selects a first element of the image; the controller 1 50A selects a lookup 

4 table 130 appropriate to that image element. It is assumed that the 

5 imaging system has already been characterized as previously described. 

6 Thus each lookup table 1 30 already contains coordinates.of locations at 

7 which the image is to be sampled before the offsets required by 

8 convolution are applied. The samples are taken at unequal intervals, the 

9 intervals being selected to effect an implicit weighting of the samples 

10 corresponding to a predetermined weighting function. In practice the 

1 1 lookup tables .1 30 are typically implemented as read-only memories 

12 (ROMs). The controller 1 50A then steps through the selected lookup table 

13 130. The arithmetic; logic unit 120A sums the coordinates from the: lookup 

14 table with the spatial offsets to generate a pixel Goordinate.,/The,pixel 

15 coordinate is used to retrieve a pixel value. PixeLvalu.es are accumulated 

16 by 1 20B. The accumulated result is passed to 1 20C which multiplies the^ : 

17 value by a scale factor from 1 40, selected to correspond to the selected 1 

18 lookup table 1 30. The output of 1 20C is thus an approximation to a 

19 weighted integral for that image element which comprises a correction for 

20 that pixel, i The weighted integral is subtracted from the pixel and the 

21 result is stored in the image memory 1 1 0B. The controller 1 50C selects a 

22 next image element and the above steps are repeated to compute 

23 weighted integrals associated with each of the plurality of image elements. 

24 As described above, each element of the corrected image may be 

25 computed as soon as the weighted integral for that image element has 

26 been found or the entire corrected image may be computed after all 

27 weighted integrals for all the image elements have been found. 

28 The apparatus 1 00 may be implemented as a so-called application 

29 specific integrated circuit (ASIC) using any commercially available "ASIC 

30 toolbox" or digital signal processor (DSP). The integrated circuit 

31 implementation is preferred because permits the apparatus 100 to be 

32 physically housed within a digital camera. 

33 Corrected pixels may be stored in a separate partition of the 

34 memory 1 1 0 as shown in Figure 4, or corrected pixels may be returned to 
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1 the input (i.e., ufrcorrected ) image "memory location "ereby doing the 

2 correction "in place." This savings in image memory obtains because 

3 each pixel is corrected by only a fraction of its full-scale dynamic range. 

4 Thus, corrected and uncorrected images are approximately the same for 

5 the purpose of computing the correction. The two-memory and one- 

6 ; memory techniques are well known as point-simultaneous and point- 

7 successive iterative methods, respectively, to those who solve dominant- 

8 diagonal sets of simultaneous linear equations. See Jansson, previously 

9 cited, pp. 85-89.- Besides saving memory, the modified method converges 

10 faster when iterations are required. Faster convergence results in less 

1 1 error, even when only one iteration is employed, as in Equation (9). 

12 However, point^successive methods are known to introduce solution 

13 ^ assymetry when only partially converged because processing is 

14 sequential. In the present application this assymetry is readily mitigated 

15 by processing pixels in the order given by address^bit reversal: See 

16 Jahsson, previously cited, pp. 85-89. 

17 Example 1 In order to correct the image for stray-rlight 

18 contamination, it is necessary to compute an object estimate according to 

19 Equation(9). This equation contains the matrix product s^i, which, 

20 according to the foregoing, may also be expressed with matrix elements 

21 proportional to the weighted sum of image pixels i(x',y), * x 

22 ZZ^U-^-/)^/), (15) ' 

* / 

23 where the sum is taken over all the image coordinate values x' and y' for 

24 each image pixel to be corrected. 

25 According to the method of the present invention, for each such 

26 pixel, this sum may be replaced by the unweighted sum- 

27 2>V,/). (16) 

selected set 
of (*',/) 
" coordinates 

28 This sum is riot taken over all the image coordinate values, but rather 

29 merely over a certain set of image pixel coordinates chosen according to 

30 the algorithm shown in Figure 2, or by some other appropriate means 

31 satisfying the requirements of implicit weighting as previously described. 
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1 In keeping wiWie coordinate requirements for a cnHrete convolution, 

2 each pixel to be corrected will require different x and y offsets between the 

3 coordinates of i(x,y) and the sampling pattern corresponding to s^x.y). 

4 This may be clarified by reference to the example illustrated in Fig. 1A. 

5 Figure 5 depicts a highly simplified case of selected ordinate 

6 sampling of the image. A 12-pixel-high by 16-pixel-wide image is shown 

7 as a checkerboard-like array of pixel squares. A hypothetical sample set 

8 of selectedrordinate coordinates are plotted by sample number (1 through 

9 14) in the pixels selected to effect the desired implicit weighting. Note that 

10 sample number 14 actually lies outside the 1 2 by 1 6 image. The sample- 

1 1 set coordinates are listed in the table to the right.. For purposes of this 

12 example, it is assumed that the sampling pattern to effect this weighting 

13 could have been chosen by means of the algorithm of Figure 2 or by any 

14 other appropriate means? Plot coordinates are given relative to both pixel ! 
15 :. . coordinates of the image /, and to the coordinates of. the spread function ; 

16 s b . Image values at the selected locations are then summed to effect a 

17 - correction to the pixel identified by the circle, which corresponds to s^0,0). 

1 8 The image depicted in Figures 5 and 6 is to be understood as comprising 

19 12 isoplanatic 4-pixel by 4-pixel patches identified in Figures 5 and 6 by 

20 the heavier outlines. All pixels within the 4 x 4 square containing, the circle 

21 will be corrected by employing a sum of samples selected by using the 

22 same pattern. That is, the pixels selected will bear the same,spatial 

23 relationship to each other and to the pixel to be corrected, for each of the 

24 image pixels to be thus corrected within the shaded isoplanatic patch. 

25 To f urther clarify the sample shifting necessary in order to compute 

26 the convolution, Figure 6 shows the image samples selected to correct an 

27 adjacent pixel in the same isoplanatic patch. It is noteworthy that this shift 

28 has moved sample 1 1 outside the bounds of the available image 

29 information, as is sample 14 both in the present and previous cases. A 

30 question thus arises as to the value assigned to the missing pixels. As 

3 1 has been explained in a previous section, one expedient is to employ the 

32 value of the closest available pixel, thereby implicitly assuming that the 

33 image would be slowly varying across its boundary had the missing pixels 

34 also been available. 
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1 It shouldW noted that an identical procedure TWbllowed for all 

2 image pixels, even those outside the shaded isoplanatic patch. In the 

3 more general case of a shift-variant spread function, however, twelve 

4 different sample patterns would be required, a different sample pattern to 

5 be used for each isoplanatic patch. In the general case of the example 

6 illustrated, the pattern shown in Figures 5 and 6 would only be useful for 

7 , the shaded patch. 

8 Applications to Volume Reconstruction As previously stated, the 

9 present method may be applied to three-dimensional situations. One or 

10 more image target volume element values of a volume reconstruction may 

1 1 contain unwanted contributions, such as might occur as a result of 

12 scattering phenomena. Examples of methods that yield such 

13 reconstructions include x^ray computed tomography (GT), magnetic- 

14 resonance imaging (MRI), and positron emission tomography (PET). In 

15 **»» thesev reconstructions, target volume element values might-represent x-ray 

16 ° absorbance, x-ray density, absorption coefficient or other physical 

17 property, or the values mightrepresent flux emitted from object volume 

18 portions corresponding to the image target volume elements. In these 

19 cases the target volume data may be processed by three-dimensional 

20 selectedrordinate convolution in a manner fully analogous to the £wo- 

21 dimensional selected-ordinate convolution in order to remove the 

22 unwanted contributions (e.g., stray or misdirected flux), or to effect 

23 enhancement of-some image-volume property such as .resolution. 

24 Under these circumstances, it should be appreciated that the two 

25 dimensional selected-ordinate method of computing a discrete 

26 convolution, and its extension to the shift-variant case via convolution 

27 within isoplanatic patches, may be readily extended to an arbitrary integral 

28 number of dimensions merely by changing the dimensionality of the 

29 associated integrals and sums in a straightforward way. In a three- 

30 dimensional case of tomography, all the volume elements, or voxels; of 

31 the object have potential to contribute to each single target voxel of an 

32 image target volume via weighted sums. An approximation to each voxel 

33 of the target image volume may be obtained by extending the selected- 

34 ordinate method to three dimensions. Thus, the selected-ordinate 
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1 samples in thwiiree-dimensional application are a"ecific set of object 

2 voxel values. The sum of these samples implicitly contains the weighting 

3 function required to compute a single target voxel value by virtue of the 

4 volume density of those selected-ordinate samples, in a manner fully 

5 analogous to the two-dimensional case. In accordance with this concept, 

6 selected-ordinate convolution may be applied to the image target volume 

7 to remove unwanted contributions from each of its voxel values. 

8 Besides the removal of unwanted contributions such as stray flux, 

9 the present method of selected-ordinate convolution may be employed to 

10 filter volume data, such as those data resulting from computed 

11 tomography or some other volume reconstruction method. Such filtering 

12 may enhance sharpness, resolution, or for reveal some .other desired 

13 property of an objects In these filtering applications, sampiingipattems 

14 corresponding to convolution weightings appropriate to the task are 
15. , employed.. 

16 Analogous to the two-dimensional case, the form or magnitude of 

17 the weighting function used for three dimensions may vary slowly over the 

18 image volume such that a condition of three-dimensional shift variance 

19 applies. A uniform region, in this case a uniform volume, may be defined 

20 that is analogous to the two-dimensional isoplanatic patch. In suph a 

21 uniform volume the form or magnitude of the weighting function varies less 

22 than some predetermined amount. Within that uniform region, selected- 

23 ordinate convolution is deemed an adequate approximation to shift-variant 

24 processing. In order to process the entire image volume, the region (in 

25 this case, the volume) is subdivided into a plurality of such uniform 

26 regions. Each uniform region will have an associated set of weights and a 

27 corresponding sampling pattern. Again, the processing is analogous to 

28 that employed in a two-dimensional situation, so that selected-ordinate 

29 convolution may be carried out for each uniform region or volume by using 

30 a sample pattern identified with that volume, thereby to carry out a shift- 

3 1 variant processing over the entire image volume. 

32 Example 2 Higher Dimensions: An Example from Four 

33 Dimensions Applications in four and higher dimensions are also practical. 

34 It is anticipated that the processing advantage (i.e., reduced computational 
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1 requirement) ofWe selected-ordinate method over cowentional methods 

2 would increase with increasing dimensionality. To cite an example of a 

3 four-dimensional case, consider the time variation of a series of 3-D 

4 tomographic reconstructions, in this example, time comprises the fourth 

5 dimension. Each successive reconstruction in the series could contain 

6 contributions from the previous reconstructions. Thus, four-dimensional 

7 filteringiby the selected- ordinate method would be required to process the 

8 * four-dimensional data to remove unwanted contributions originating from 

9 various locations, both spatial and temporal. In analogy to the two- and 

10 three-dimensional cases, shift-variant processing may be accomplished by 

1 1 dividing the* image into uniform regions having smaller dimensions, but the 

12 same 4 dimensionality as the image. 

13 In the context of image correction, the term "stray:flux"<or 

14 "smearing" is used to describe corruption effects upon an image; It is to 

15 be understood that the term "smearing" includes effects described by the 

16 term "stray flux" and also includes "temporal smearing". 

17 ;ln the context of image formation, the term "misdirected" is used to 

18 describe flux that is analogous to the "stray flux" or "smearing". The term 

19 "misdirected" also includes "temporally misdirected". 

20 Example 3 > Application to Finite Impulse Response Filtering (FIR) 

21 in One Dimension It should also be appreciated that a orte-dimensional 

22 selected-ordinate filtering method is also useful. The previously discussed 

23 spectrometry determination of single, separate colorimetric X, Y, and Z 

24 values employed selected-ordinate sampling but not the selected-ordinate 

25 convolution of the present method. Application of discrete convolution, 

26 and its analog, by selected-ordinate convolution in one dimension finds 

27 specific application in circuits to carry out finite-impulse-response filtering 

28 of electronic signals. Such filtering could include cases where the filter 

29 transfer function varies with time: a shift variant case. High-speed tapped 

30 analog delay line finite-impulse-response filters may be designed in 

31 accordance with the selected-ordinate method. By reducing the number of 

32 samples and eliminating multiplication operations, cost and complexity are 

33 reduced while speed of filter response is increased. 
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1 ExampW Mixed Dimensions Selected-ordMte convolution also 

2 has practical application to image formation, image filtering and image 

3 correction in cases where the object and image differ in their 

4 dimensionality- 1 Suppose, for example, a two-dimensional target plane 

5 receives contributions from a number of volume element portions of .an 

6 object volume. Telocentric microscopy of a volume of fluorescing light 

7 sources is such an example. Each image target element receives flux 

8 primarily from emitters along a line parallel to the. optical axis. It also 

9 receives flux from adjacent lines owing to light scattering in the object, to 

10 finite depth of focus, to optical, diffraction, and other phenomena 

1 1 characterized by a volume point-spread function. Also, some of the flux 

12 may be absorbed in the object. Thus the flux arriving at a given image- 

13 plane target element is affected by. a plurality of object volume portions. 

14 To an approximation, the equation that describes this situation is 

15 i(x,y)= jls(Ky,x'y,^Xy,*)dx'd/dz'. *r (17) 

16 In this equation, o(x',y',z) is the flux emitted byWie" object volume, i{x,y) is 

17 ' the projected image, and s(x,y;x',y',z) is the projection iterator that 

18 describes the physics of this scenario. The quantities x and y are spatial 

19 cartesian coordinates of the image, and x',y', and z 1 are Volume Y 

20 coordinates within the target volume. As before, if the function 

21 s{x,y;x',y',z) cannot be reformulated as a convolution over-the entire 

22 object, it may be possible to do so over a smaller uniform volume. When 

23 this situation applies, as should almost always be the case, selected- 

24 ordinate convolution can be Used to compute i{x,y) values from o(x',y',z) 

25 values by applying' Equation (1 7) above. Accordingly samples of o{x',y',z) 

26 are chosen according to the selected-ordinate method. Accordingly, the 

27 i{x,y) values may be processed via selected ordinates to remove 

28 unwanted contributions from undesired portions of o(x',y',z). Also, image 

29 i(x,y) may be filtered via selected-ordinate convolution to enhance 

30 resolution or otherwise improve or alter its properties in a desired fashion. 

31 In some cases, it will be desired to selectively sample computed or 
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the object elements in accordancHvil 



1 assumed valuePor the object elements in accordancWvith the selected- 

2 ordinate method in order to effect the correction. 

3 Analogously, selected-ordinate convolution may be applied to 

4 objects and corresponding images having different dimensionalities. 

5 Applications include obtaining images from objects^ filtering, correcting, 

6 and otherwise' processing images having arbitrary dimensionality. 

7 It is to be understood that in when applying the selected-ordinate 

8 method in the case of mixed dimensionality, discrete summations take the 

9 place of the integral of Equation (17) in the foregoing, and that the role of 

10 the weighting described by s(x,y;x',y ',z) is assumed by the spacing of the 

11 volume element samples. As in the two-dimensional casey*the weighting 

12 is implicit in- the sample spacing: It is to be anticipated that special 

13 purpose processors, similar to those employed in the two-dimensional . 

14 case rcan also be applied to effect selected-ordinate convolution for 

15 applications having arbitrary and mixed dimensionalities. ; 

' — «.i 

16 Iteration It should be appreciated that whereas therp/eferred mode 

17 of image cor r£ctioifcky«,select.ed-ordinate processing employs only one [ 

18 cycle of.an^therwise iterative computation exemplified by Equation (9), vi 

19 iteration may yield a-more highly corrected image. Accordingly, an . 

20 iterative, method,, as described by Equation (4), may be applied to images 

21 having any dimensionality. In such applications, it may also be desired to 

22 alter the normalized value of s during some of or all of the cycles of 

23 iteration such that s P does not equal I (the identity matrix) for those 

24 iterations. As is typical for iterative methods, cycles of iteration may 

25 proceed until a fixed amount of time has elapsed, a fixed number of cycles 

26 have been executed,, or the correction-term magnitude becomes smaller 

27 than some predetermined stopping value. Other criteria may also be 

28 used, either alone or in combination with one or more of the aforesaid 

29 criteria. 
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What is claimed is: 
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1 . A method of correcting an image of an object, the image 



5 comprising one or more dimensions and having a plurality of image 

6 elements and the image being corrupted by smearing, to remove the effect 

7 of the smearing on the image, the method comprising the steps of: 
8 

9 a) sampling portions of the image of the object at unequal intervals , 

10 the intervals being selected to effect an implicit weighting of the samples 

1 1 corresponding to a predetermined weighting function; 

12 .••.■:-> - 

13 b) summing the implicitly weighted samples to create an 

14 approximation to a weighted integral; 
15 

16 c) repeating steps a) - b) to compute weighted integrals associated 

17 with each of the plurality of image elements thereby to obtain a correction 

18 to the image. 
19 

20 2. The method of claim 1 further comprising, before step a), the > 

21 step of determining a weighting function that characterizes the^smearing. 

22 

23 3. The method of claim 1 further comprising: 

24 d) applying the correction to the image to remove therefrom the 

25 smearing corruption to create a corrected image. 

26 

27 4. The method of claim 3 further comprising the steps of: 

28 e) applying an additional correction term containing a sum of 

29 implicitly weighted samples to each element of the corrected image to 

30 update the corrected image; 

31 f) repeating step e) until a predetermined stopping criterion has 

32 been satisfied. 

33 
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1 5. The mwiod of Claim 4, where the image isTBed as a first 

2 approximation to the corrected image. 

3 

4 6. The method of Claim 4, wherein the number of dimensions of the 

5 image differs from the number of dimensions of the object. 

6 • •• • ., 

7 7. A method of obtaining values of image elements, the image 

8 having a plurality of dimensions, in an flux imaging system characterized 

9 by a portion of the image flux being misdirected, the method comprising 

10 the steps of: 

11 a) sampling portions of an object at unequal intervals, the intervals 

12 being selected to effect an implicit weighting of the samples corresponding 

13 to a predetermined weighting f unction; 

14 b) summing the implicitly weighted samples to create an 

15 approximation to a weighted integral; 

16 ■ c) repeating steps a) - b) to compute weighted integrals associated 

17 with a plurality of image elements thereby to obtain the values of the 

18 image elements. 
19 

20 8. The method of claim 7 further comprising, before step a), 

21 determining a weighting function that characterizes the misdirection. 

22 

23 9. The method of Claim 7, wherein the number of dimensions of 

24 the image differs from the number of dimensions of the object. 

25 

26 1 0. A method of correcting an image of an object, the image 

27 having a plurality of image elements and the image being corrupted by 

28 stray flux, to remove the effect of stray flux on the image, the method 

29 comprising the steps of: 

30 a) sampling portions of the image of the object at unequal intervals, 

31 the intervals being selected to effect an implicit weighting of the samples 

32 corresponding to a predetermined weighting function; 

33 b) summing the implicitly weighted samples to create an 

34 approximation to a weighted integral; 

30 



1 c) rep Jfcig steps a) - b) to compute weigh^hntegrals associated 

2 with each of the plurality of image elements thereby to obtain a correction 

3 to the image. 

4 • ■ •• • . 

5 11. The method of claim 10 wherein the image comprises a 

6 plurality of uniform regions, wherein the weighting function within a uniform 

7 region does not vary by more than a predetermined amount. 

8 ' . 

9 12. The method of claim .1 1 wherein each uniform region of the 

10 object or its corresponding image is associatedwith a different sampling 

11 pattern. : 

12 

13 13. The method of claim 1 0 further comprising: 

14 d) applying the correction to the image to remove therefrom the 

15 stray-flux corruption. 

17 14. A method of filtering an image , the method comprising the 

18 steps of: 

19 a) sampling portions of the image at unequal intervals, the intervals 

20 being selected to effect an implicit weighting of the samples corresponding 

21 to a predeterminedtweighting function; 

22 b) summing the implicitly weighted samples to create an 

23 approximation to a weighted integral; 

24 c) repeating steps a) - b) to compute weighted integrals associated 

25 with a plurality of image elements thereby to compute a filtered image. 

26 . . • . 

27 1 5. The method of claim 14 wherein the image comprises a 

28 plurality of uniform regions, wherein the weighting function within a uniform 

29 region does not vary by more than a predetermined amount. 
30 

31 16. The method of claim 15 wherein each uniform region of the 

32 object or its corresponding image is associated with a different sampling 

33 pattern. 
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n^^od of obtaining values of image ellfl^ 



1 17. A melnod of obtaining values of image eleTfrents, the method 

2 r comprising the steps of: 

3 

4 a) sampling portions of an object at unequal intervals, the intervals 

5 being selected to effect an implicit weighting of the samples corresponding 

6 to a predetermined weighting function; 

7 b) summing the implicitly weighted samples to create an 

8 approximation to a weighted integral; 

9 c) repeating steps a) - b) to compute weighted integrals associated 

10 with a plurality of image elements thereby to obtain the values of the 

11 image elements. 
12 

13 18. The method of claim 17, wherein thehumber of dimensions of 

14 the image differs from the number of dimensions of the object; 
15 

16 1 9. The method of claim 17 wherein the image comprises a 

17 plurality of uniform regions, wherein the weighting function within a uniform 

18 region does not vary by more than a predetermined amount. . 

19 -<■.•: . • .- • • ■ * v- 

20 20. The method of claim 19 wherein each uniform region:of the 

21 object or its corresponding image is associated with a different sampling 

22 pattern. >■ 
23 

24 21. A method for imaging an object having a predetermined 

25 number of discrete portions thereon, the object having flux emanating 

26 therefrom, by generating a signal representative of the total incident flux 

27 on each area element of a target plane, each area element corresponding 

28 to a portion of the object, the signal being representative of the sum of flux 

29 emanating from the corresponding portion of the.object plus flux scattered 

30 from scattering volumes external to the object, the method comprising the 

31 steps of: 

32 a) characterizing the flux scattered by the scattering volumes from 

33 each portion of the object to create a weighting function representing the 

34 total flux incident on a particular target-plane area element; 

32 



1 b) sam^fcg the portions of the object at une^fcl intervals, the 

2 intervals being selected to effect an implicit weighting of the samples 

3 corresponding to the weighting function created in step a); 

4 c) summing the implicitly weighted samples to create an 

5 approximation to a weighted integral; 

6 d) repeating steps a) - c) for each target-plane element to compute 

7 weighted integrals associated with a plurality of target-plane elements 

8 thereby to compute an image correction; 

9 e) employing the image correction thus obtained to remove from 

10 target-plane flux signals the contribution to such signals caused by 

11 scattered flux. 

12 - 

13 22. The method of claim 21 , wherein the target plane comprises a < 

14 plurality of smaller isoplanatic patches,- and wherein such patches may be 

15 associated with different weighting functions. 
16 

17 23. The method of claim 22 wherein each isoplanatic patch has a 

18 sample pattern associated therewith. 

19: ~ 

20 • . • „ 24. The method of claim 22 wherein the size of each isoplanatic 

21 patch is associated with the spatial rate of change of the weighting 

22 function. 

23 

24 25. The method of claim 21 wherein the weighting function within 

25 each isoplanatic patch does not vary by more than a predetermined 

26 amount. 

27 

28 26. The method of claim 21 , wherein the sample intervals of step b) 

29 for sampling the portions of the object at unequal intervals are obtained 

30 from an image pyramid. 
31 

32 27. An apparatus for correcting an image of an object, the image 

33 having a plurality of image elements and the image being corrupted by 
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1 stray flux, to refl^ve the effect of stray flux on the irn^e, the apparatus 

2 comprising: 

3 a) an image memory having one or more partitions for storing an 

4 input image and a corrected image; 

5 b) an arithmetic logic unit for performing arithmetic computations; 

6 c) a plurality of look up tables for storing selected-ordinate image 

7 coordinate values; 

8 d) a controller for addressing the image memory and selecting 

9 lookup tables; 

10 where, in operation 

11 1) an image of the object is stored in the image memory; : 

12 2) the controller selects an element of the image; 

13 3) the controller selects a lookup table, and a corresponding scale 

14 factor, appropriate to that image element, each lookup table containing 

15 information to generate a pattern of sample coordinates of .the image of 

16 the object, the samples being at unequal intervals, the intervals being 

17 selected to effect an implicit weighting of the samples corresponding to a 

18 predetermined weighting function; 

19 4) the controller steps through the lookup table and the arithmetic 

20 logic unit sums the implicitly weighted samples and multiplies thejn by the 

21 scaling factor to create an approximation to a weighted integral for that 

22 image element; 

23 5) the weighted integral is stored in the image memory; 

24 6) the controller selects a next element and steps 3) - 5) are 

25 repeated to compute weighted integrals associated with each of the 

26 plurality of image elements, 

27 thereby to obtain a corrected image. 
28 

29 28. The apparatus of claim 27 wherein 

30 the weighted integral created in step 5) is subtracted from the 

31 image value and the resulting value is stored in the image memory. 

32 

33 29. A computer-readable medium containing instructions for 

34 providing program control for a computer controller to remove the effects 

34 



of the spatial Jfcaring on an image of an object, th^nage comprising 
one or more dimensions and having a plurality of image elements, by 
implementing the steps of: 

a) sampling portions of the image of the object at unequal intervals , 
the intervals being selected to effect an implicit weighting of the samples 
corresponding to a predetermined weighting function; 

b) summing the implicitly weighted samples to create an 
approximation to a weighted integral; 

c) repeating steps a) - b) to compute weighted integrals associated with each 
of the plurality of image elements thereby to obtain a correction to the image. 
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3 Method, Program and Apparatus for Efficiently Removing 

4 Stray-Flux Effects by Selected-Ordinate Image Processing 

5 Abstract of the Invention 

6 The present invention is directed toward correctinga corrupted 



7 image of an object to remove the effects of smearing or misdirected flux in 

8 the image, wherein the image comprises one or more dimensions and has 

9 a plurality, of image elements. The method of the present invention 

10 comprises the steps of: 

11 a) sampling portions of the image of the object at unequal intervals, 

12 the intervals being selected to effect an implicit weighting of the samples 

13 corresponding to a predetermined weighting function; 

14 b) summing the implicitly weighted samples to create an 

15 approximation to a weighted integral; 

16 c) repeating steps a) - b) to compute weighted integrals associated 

17 with each of the plurality of image elements thereby to obtain a correction 

18 to the image. 

19 The present invention is also implemented as a computer readable 

20 medium for controlling a computer to implement the steps of the method 

21 and in apparatus form. 
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ALGORITHM TO GENERATE SELECTED-ORDINATE 
COORDINATES 



DO FOR ALL ISOPLANATIC PATCHES 



DO FOR A REPRESENTATIVE PIXEL IN EACH PATCH 



PLACE PIXEL COORDINATE LOCATION OF WEIGHT- 
FUNCTION MAXIMUM ON SELECTED-ORDINATE LIST 

DO FOR MONQTONIGALLY INCREASING REPRESENATATIVE- 
PATCH-ASSOCIATED WEIGHT FRACTIONS UNTIL FULL 
WEIGHT IS REALIZED 



DO FOR ALL IMAGE PIXELS 



COMPARE LOCAL SAMPLE DENSITY WITH FRACTION 
OF REPRESENTATIVE-PATCH ASSOCIATED WEIGHT 

IF LOCAL SAMPLE DENSITY IS INADEQUATE & NO 
SAMPLE HAS BEEN PLACED AT THIS LOCATION 



ADD THIS PIXEL'S COORDINATE LOCATION TO 
SELECTED- 
ORDINATE LIST; FOR THIS PATCH 



ALL ORDINATES FOR THIS PATCH HAVE BEEN SELECTED 

DETERMINE SCALE FACTOR NEEDED TO RELATE SELECTED 
IMAGE PIXEL SUM TO STRAY-LIGHT ESTIMATE 



Figure 2 



STRAY-LIGHT CORRECTION ALGORITHM 

DO FOR ALL ISOPLANATIC PATCHES 
DO FOR ALL PIXELS IN EACH PATCH 

CLEAR STRAY-LIGHT ACCUMULATOR 

DO FOR ALL PATCH SELECTED-ORDINATE 
COORDINATES 



ADD PATCH-SPECIFIC OFFSET 
TO SELECTED-ORDINATE COORDINATE 

ADD TO ACCUMULATOR THE PIXEL VALUE 
|_ AT THE OFFSET-CORRECTED COORDINATE 

SCALE ACCUMULATED STRAY-LIGHT SUM 

SUBTRACT STRAY-LIGHT ESTIMATE FROM IMAGE 
PIXEL 



Figure3 




Figure 
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ELECTED-ORDINATE SAMPLE 
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i Coordinate 
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